rm(list = ls())
source("C:/Users/mbofi/Dropbox/CeMSIIS/GitHub/Allocation/case-study/aux_functions.R") #local
# setwd("C:/Users/mbofi/Dropbox/CeMSIIS/GitHub/Allocation/case-study")#local
# setwd("~/GitHub/Allocation/case-study")
library(tidyverse)
#> Warning: package 'tidyverse' was built under R version 4.2.1
#> ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
#> ✔ ggplot2 3.4.1 ✔ purrr 0.3.4
#> ✔ tibble 3.1.8 ✔ dplyr 1.0.10
#> ✔ tidyr 1.2.1 ✔ stringr 1.4.1
#> ✔ readr 2.1.2 ✔ forcats 0.5.2
#> Warning: package 'ggplot2' was built under R version 4.2.3
#> Warning: package 'tibble' was built under R version 4.2.1
#> Warning: package 'tidyr' was built under R version 4.2.1
#> Warning: package 'readr' was built under R version 4.2.1
#> Warning: package 'dplyr' was built under R version 4.2.1
#> Warning: package 'stringr' was built under R version 4.2.1
#> Warning: package 'forcats' was built under R version 4.2.1
#> ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
#> ✖ dplyr::filter() masks stats::filter()
#> ✖ dplyr::lag() masks stats::lag()
library(dplyr)
library(plyr)
#> Warning: package 'plyr' was built under R version 4.2.1
#> ------------------------------------------------------------------------------
#> You have loaded plyr after dplyr - this is likely to cause problems.
#> If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
#> library(plyr); library(dplyr)
#> ------------------------------------------------------------------------------
#>
#> Attaching package: 'plyr'
#>
#> The following objects are masked from 'package:dplyr':
#>
#> arrange, count, desc, failwith, id, mutate, rename, summarise,
#> summarize
#>
#> The following object is masked from 'package:purrr':
#>
#> compact
library(NCC)
#> Registered S3 methods overwritten by 'registry':
#> method from
#> print.registry_field proxy
#> print.registry_entry proxy
#> Warning: package 'memoise' was built under R version 4.2.1
library(mmtable2)
#>
#> Attaching package: 'mmtable2'
#>
#> The following object is masked from 'package:tidyr':
#>
#> table1
library(gt)
#> Warning: package 'gt' was built under R version 4.2.1
# OPTIMAL ALLOCATION USING OLD NOTATION
##########################################
# Simulation function
# modification of sim_designs() function
# also returning the matrix (r_is)
##########################################
sim_designs_rr <- function(r1,r2,mu0,mu1,mu2,N,alloc="sqrt",trend="stepwise",sl=0.2){
r3 = 1-r1-r2
if(r1 == 1){
if(alloc == "one"){
r22 <- r11 <- r01 <- r1/3
r02 <- r12 <- r23 <- r03 <- 0
}
if(alloc != "one"){
r01 <- sqrt(2)/(2+sqrt(2))
r22 <- r11 <- (1-r01)/2
r02 <- r12 <- r23 <- r03 <- 0
}
}else{
if(alloc == "opt"){
r11 <- r01 <- r1/2
r2_opt <- f(r1=r1,r2=r2)
r02 <- r2_opt[3]
r12 <- r2_opt[4]
r22 <- r2_opt[5]
r23 <- r03 <- r3/2
}
if(alloc == "one"){
r11 <- r01 <- r1/2
r22 <- r12 <- r02 <- r2/3
r23 <- r03 <- r3/2
}
if(alloc == "sqrt"){
r11 <- r01 <- r1/2
r02 <- r2*sqrt(2)/(2+sqrt(2))
r22 <- r12 <- (r2-r02)/2
r23 <- r03 <- r3/2
}
}
# c(r11,r01)
n11 = round(r11*N)
n01 = round(r01*N)
n22 = round(r22*N)
n12 = round(r12*N)
n02 = round(r02*N)
n23 = round(r23*N)
n03 = round(r03*N)
# c(r11,r01,r22,r12,r02,r23,r03)
# c(n11,n01,n22,n12,n02,n23,n03)
means = c(mu0,mu1,mu2)
treatment = c(
sample(c(rep(1,n11),rep(0,n01))),
sample(c(rep(2,n22),rep(1,n12),rep(0,n02))),
sample(c(rep(2,n23),rep(0,n03)))
)
Nsim = length(treatment)
if(r1==1){
treatment = sample(treatment)
period = rep(1,Nsim)
}else{
period = c(
rep(1,n11+n01),
rep(2,n22+n12+n02),
rep(3,n23+n03)
)
}
if(trend=="stepwise"){
response = rnorm(Nsim,
mean=means[treatment[1:Nsim]+1]+sw_trend(cj=period[1:Nsim], lambda=sl),
sd=1)
}
if(trend=="linear"){
response = rnorm(Nsim,
mean=means[treatment[1:Nsim]+1]+linear_trend(j=1:Nsim, lambda=sl, sample_size=c(0,Nsim)),
sd=1)
}
data = data.frame(response,treatment,period)
if(r1==1){
ss = matrix(c(n22,n11,n01,0,n12,n02,n23,0,n03), nrow=3)
}else{
ss = matrix(c(0,n11,n01,n22,n12,n02,n23,0,n03), nrow=3)
}
r_matrix = matrix(c(0,r11,r01,r22,r12,r02,r23,0,r03), nrow=3)
r_periods = c(r1,r2,r3)
return(list(data=data,ss=ss,r_matrix=r_matrix,r_periods=r_periods))
}
# OPTIMAL ALLOCATION USING THE NEW NOTATION
##########################################
eq_p22 <- function(p22,r1=0.1,r2=0.8){(Power(-1 + p22,3)*(-1 + 2*r1))/((-1 + 2*p22)*(-2 + p22*(7 + p22*(-15 + p22*(19 + 2*p22*(-7 + 2*p22))))))-r2}
f_p=Vectorize(function(r1,r2) {
p22=uniroot(eq_p22,c(0,r2),r1=r1,r2=r2)$root
r22 = r2*p22
# p02 <- 1/(2-2*p22)
p02 <- 1/(2-2*p22)-p22
r02 <- p02*r2
p12 <- 1-p02-p22
r12 <- p12*r2
# r12 <- r2- r02- r22
sol=c(r1,r2,r22,r12,r02)
sol
})
##########################################
# Case study
##########################################
# original study - obtained means
mean_control = 17.3/3.5
mean_arm1 = 66.2/3.5
mean_arm2 = 72.3/3.5
##########################################
# design 3: three-period design (symmetric design)
##########################################
N = 92
N1 = round(N/3)
N3 = round(N/3)
N2 = N-N1-N3
c(N1,N2,N-N1-N2)
#> [1] 31 30 31
alloc_str <- c("one","opt","sqrt")
Optimal allocation using the old parametrization
m <- sim_designs_rr(r1=N1/N,r2=N2/N,
mu0=mean_control,mu1=6,mu2=6,
N=N,alloc=alloc_str[2],sl=0)
rownames(m$r_matrix) <- c("Arm 2", "Arm 1", "Control")
knitr::kable(m$r_matrix, format = "html", caption = paste(alloc_str[2]),
col.names = c("Period 1", "Period 2", "Period 3"),
row.names=T,
digits=3)
| Period 1 | Period 2 | Period 3 | |
|---|---|---|---|
| Arm 2 | 0.000 | 0.096 | 0.168 |
| Arm 1 | 0.168 | 0.096 | 0.000 |
| Control | 0.168 | 0.135 | 0.168 |
knitr::kable(t(m$r_periods),col.names = c("Period 1", "Period 2", "Period 3"))
| Period 1 | Period 2 | Period 3 |
|---|---|---|
| 0.3369565 | 0.326087 | 0.3369565 |
Optimal allocation using the new parametrization
r_solutions <- f_p(r1=m$r_periods[1], r2=m$r_periods[2])
# rownames(r_solutions[3:5]) <- c("Arm 2", "Arm 1", "Control")
knitr::kable(r_solutions[3:5],col.names = c("Period 2"),row.names=T,digits=3)
| Period 2 | |
|---|---|
| 1 | 0.096 |
| 2 | 0.096 |
| 3 | 0.135 |
knitr::kable(t(c(r_solutions[1:2], 1-sum(r_solutions[1:2]))),col.names = c("Period 1", "Period 2", "Period 3"),digits=3)
| Period 1 | Period 2 | Period 3 |
|---|---|---|
| 0.337 | 0.326 | 0.337 |
##########################################
# design 3: three-period design (non-symmetric design)
##########################################
N = 92
N1 = round(N/3)
N2 = round(2*(N-N1)/3)
c(N1,N2,N-N1-N2)
#> [1] 31 41 20
Optimal allocation using the old parametrization
m <- sim_designs_rr(r1=N1/N,r2=N2/N,
mu0=mean_control,mu1=6,mu2=6,
N=N,alloc=alloc_str[2],sl=0)
rownames(m$r_matrix) <- c("Arm 2", "Arm 1", "Control")
knitr::kable(m$r_matrix, format = "html", caption = paste(alloc_str[2]),
col.names = c("Period 1", "Period 2", "Period 3"),
row.names=T,
digits=3)
| Period 1 | Period 2 | Period 3 | |
|---|---|---|---|
| Arm 2 | 0.000 | 0.169 | 0.109 |
| Arm 1 | 0.168 | 0.087 | 0.000 |
| Control | 0.168 | 0.190 | 0.109 |
knitr::kable(t(m$r_periods),col.names = c("Period 1", "Period 2", "Period 3"),digits=3)
| Period 1 | Period 2 | Period 3 |
|---|---|---|
| 0.337 | 0.446 | 0.217 |
Optimal allocation using the new parametrization
r_solutions <- f_p(r1=m$r_periods[1], r2=m$r_periods[2])
# rownames(r_solutions[3:5]) <- c("Arm 2", "Arm 1", "Control")
knitr::kable(r_solutions[3:5],col.names = c("Period 2"),row.names=T,digits=3)
| Period 2 | |
|---|---|
| 1 | 0.169 |
| 2 | 0.087 |
| 3 | 0.190 |
knitr::kable(t(c(r_solutions[1:2], 1-sum(r_solutions[1:2]))),col.names = c("Period 1", "Period 2", "Period 3"),digits=3)
| Period 1 | Period 2 | Period 3 |
|---|---|---|
| 0.337 | 0.446 | 0.217 |
Optimal allocation using the old parametrization
# code used in the shiny app
source("C:/Users/mbofi/Dropbox/CeMSIIS/GitHub/Allocation/optimisation/comparisons/aux_functions_ncc.R")
#> Warning: package 'nloptr' was built under R version 4.2.3
#> Warning: package 'latex2exp' was built under R version 4.2.1
#> [,1] [,2] [,3] [,4]
#> [1,] 0.2 0.5 0.376248 0.2096949
results <-optr1(r1=0.3)
colnames(results)=c("r1","r2","r12","r22")
results=as.data.frame(results)
results$r02 <- 1- results$r12-results$r22
results[c(51),]
#> r1 r2 r12 r22 r02
#> 51 0.3 0.7 0.1656498 0.4318298 0.4025203
results[c(51),3:5]*0.7
#> r12 r22 r02
#> 51 0.1159549 0.3022809 0.2817642
Optimal allocation using the old parametrization
sol_r22 <- function(r1,r2){
a <- Power(8 + 36*r1 - 108*Power(r1,2) + 27*Power(r1,3) +
6*sqrt(3)*sqrt(r1*(16 - 72*r1 + 108*Power(r1,2) - 27*Power(r1,3))),
0.3333333333333333)
b <- Power(-2 + a,2) + 6*(-4 + a)*r1 + 9*Power(r1,2)
sol <- -r1 + (sqrt((b*(-1 + r1))/a) - sqrt(-(((-1 + r1)*
(-12*sqrt(3)*Power(a,1.5)*sqrt(b*(-1 + r1))*r1 +
b*(4 + 8*a + Power(a,2) - 12*(2 + a)*r1 + 9*Power(r1,2))))/(a*b))))/
(4.*sqrt(3))
sol_p <- (12*r1 + sqrt(3)*(-sqrt((b*(-1 + r1))/a) +
sqrt(-(((-1 + r1)*(-12*sqrt(3)*Power(a,1.5)*sqrt(b*(-1 + r1))*r1 +
b*(4 + 8*a + Power(a,2) - 12*(2 + a)*r1 + 9*Power(r1,2))))/(a*b)))))/
(12.*(-1 + r1))
return(c(a,b,sol,sol_p))
# 1 - r1 + sqrt(((-1 + r1)*(9*Power(r1,2) + 6*r1*(-4 + Power(8 + 36*r1 - 108*Power(r1,2) + 27*Power(r1,3) + 6*sqrt(3)*sqrt(r1*(16 - 72*r1 + 108*Power(r1,2) - 27*Power(r1,3))),
# 0.3333333333333333)) + Power(-2 + Power(8 + 36*r1 - 108*Power(r1,2) + 27*Power(r1,3) + 6*sqrt(3)*sqrt(r1*(16 - 72*r1 + 108*Power(r1,2) - 27*Power(r1,3))),
# 0.3333333333333333),2)))/Power(8 + 36*r1 - 108*Power(r1,2) + 27*Power(r1,3) + 6*sqrt(3)*sqrt(r1*(16 - 72*r1 + 108*Power(r1,2) - 27*Power(r1,3))),0.3333333333333333))/
# (4.*sqrt(3)) - sqrt(-7.333333333333333 + 8*Power(-1 + r1,2) + (43*r1)/3. - 7*Power(r1,2) -
# ((-1 + r1)*(4 - 24*r1 + 9*Power(r1,2)))/
# (12.*Power(8 + 36*r1 - 108*Power(r1,2) + 27*Power(r1,3) + 6*sqrt(3)*sqrt(r1*(16 - 72*r1 + 108*Power(r1,2) - 27*Power(r1,3))),0.3333333333333333)) -
# ((-1 + r1)*Power(8 + 36*r1 - 108*Power(r1,2) + 27*Power(r1,3) + 6*sqrt(3)*sqrt(r1*(16 - 72*r1 + 108*Power(r1,2) - 27*Power(r1,3))),0.3333333333333333))/12. +
# (sqrt(3)*Power(-1 + r1,2)*r1*Power(8 + 36*r1 - 108*Power(r1,2) + 27*Power(r1,3) + 6*sqrt(3)*sqrt(r1*(16 - 72*r1 + 108*Power(r1,2) - 27*Power(r1,3))),0.16666666666666666))/
# sqrt((-1 + r1)*(9*Power(r1,2) + 6*r1*(-4 + Power(8 + 36*r1 - 108*Power(r1,2) + 27*Power(r1,3) + 6*sqrt(3)*sqrt(r1*(16 - 72*r1 + 108*Power(r1,2) - 27*Power(r1,3))),
# 0.3333333333333333)) + Power(-2 + Power(8 + 36*r1 - 108*Power(r1,2) + 27*Power(r1,3) + 6*sqrt(3)*sqrt(r1*(16 - 72*r1 + 108*Power(r1,2) - 27*Power(r1,3))),
# 0.3333333333333333),2))))/2.
}
sol_r12 <- function(r1,r22){
(1 - r1 - sqrt(1 - r1 - 4*r22 + 4*r1*r22 + 4*Power(r22,2)))/2.
}
sol_r22(r1=.3,r2=.7)
#> [1] 2.7275108 -0.9512085 -0.6977191 -0.9967416
# sol_r12(r1=.3,r22=sol_r22(r1=.3,r2=.7)[1])
sol_r12(r1=.3,r22=sol_r22(r1=.3,r2=.7))
#> [1] -2.0385263 -0.9712281 -0.7224809 -1.0160940
Center for Medical Statistics, Informatics and Intelligent Systems, Medical University of Vienna.
[Klassifizierung: intern]
marta.bofillroig@meduniwien.ac.at